'''
Created on Jul 9, 2009

@author: mkiyer
'''

from optparse import OptionParser
import os
from veggie.app.chipseq.hpeak import HPeak
from veggie.app.chipseq.macs import MACS

#vcap_erg_chad_1 = {'VCaP_ERG_chad_1': [('30351AAXX', ('5_1',), 'export')]}
#vcap_erg_chad_2 = {'VCaP_ERG_chad_2': [('30351AAXX', ('5_2',), 'export')]}
#vcap_ku80_chad_1 = {'VCaP_Ku80_chad_1': [('30351AAXX', ('2_1',), 'export')]}
#vcap_ku80_chad_2 = {'VCaP_Ku80_chad_2': [('30351AAXX', ('2_2',), 'export')]}
#vcap_input_chad_1 = {'VCaP_input_1': [('301YWAAXX', ('4_1',), 'export')]}
#vcap_input_chad_2 = {'VCaP_input_2': [('301YWAAXX', ('4_2',), 'export')]}
#
#vcap_panh3_1 = {'VCaP_PanH3_1': [('20E46AAXX', ('6',), 'export')]}
#vcap_panh3_2 = {'VCaP_PanH3_2': [('20E6WAAXX', ('8',), 'eland_result')]}
#vcap_panh3_all = {'VCaP_PanH3_all': [('20E46AAXX', ('6',), 'export'),
#                                     ('20E6WAAXX', ('8',), 'eland_result')]}
#
#vcap_nont_sirna_erg = {'VCaP_NonT_SiRNA_ERG': [('206KMAAXX', ('4',), 'export')]}
#vcap_erg_sirna_erg = {'VCaP_ERG_SiRNA_ERG': [('206KMAAXX', ('7',), 'export')]}
#vcap_nont_sirna_ar = {'VCaP_NonT_SiRNA_AR': [('206KMAAXX', ('6',), 'export')]}
#vcap_erg_sirna_ar = {'VCaP_ERG_SiRNA_AR': [('206KMAAXX', ('8',), 'export')]}
#
bt474_t = {'BT474_H3K4me3_E2': [('315A9AAXX', ('1',), 'export')]}
#bt474_c = {'BT474_H3K4me3_EtOH': [('315A9AAXX', ('2',), 'export')]}
mcf7_h3k4_t = {'MCF7_H3K4me1_E2': [('315A9AAXX', ('3',), 'export')]}
#mcf7_h3k4_c = {'MCF7_H3K4me1_EtOH': [('315A9AAXX', ('4',), 'export')]}
mcf7_foxa1_t = {'MCF7_FOXA1_E2': [('315A9AAXX', ('5',), 'export')]}
#mcf7_foxa1_c = {'MCF7_FOXA1_EtOH': [('315A9AAXX', ('6',), 'export')]}

vcap_foxa1 = {'VCaP_FoxA1': [('20A9RAAXX', ('2',), 'export')]}

vcap_erg = {'VCaP_ERG': [('30351AAXX', ('5_1', '5_2',), 'export')]}
vcap_ku80 = {'VCaP_Ku80': [('30351AAXX', ('2_1', '2_2',), 'export')]}
vcap_input = {'VCaP_input': [('301YWAAXX', ('4_1', '4_2',), 'export')]}
vcap_dnapk = {'VCaP_DNAPKcs': [('30351AAXX', ('3_1', '3_2'), 'export')]}
vcap_ph2ax = {'VCaP_pH2AX': [('30351AAXX', ('4_1', '4_2'), 'export')]}
vcap_control = {'VCaP_control': [('301YWAAXX', ('4_1', '4_2',), 'export'),
                               ('20E46AAXX', ('6',), 'export'),
                               ('20E6WAAXX', ('8',), 'eland_result'),
                               ('206KMAAXX', ('7',), 'export')]}
vcap_r1881_erg = {'VCaP_R1881_ERG': [('20AV3AAXX', ('8',), 'eland_result')]}

def find_peaks(peakfinder, treatmentdict, controldict=None, cwd=None):
    # run peak finder on all combinations on treatments/controls
    for treatment_name, treatments in treatmentdict.iteritems():
        if controldict:
            for control_name, controls in controldict.iteritems():
                peakfinder.setUp(treatment_name, treatments, control_name, controls, cwd=cwd)
                peakfinder()
                peakfinder.tearDown()
        else:
            peakfinder.setUp(treatment_name, treatments, cwd=cwd)
            peakfinder()
            peakfinder.tearDown()

if __name__ == '__main__':
    # setup command line parser
    parser = OptionParser()
    parser.add_option("--flowcell", dest="flowcell")
    parser.add_option("--lane", dest="lane")
    
    
    
    
    usage = "usage: %prog [options] elandfile"
    optparser = OptionParser()
    optparser.add_option("--hpeak",
                         action="store_true", dest="hpeak", default=False,
                         help="create the pseudo-BED format needed by hpeak")
    optparser.add_option("--format", action="store", dest="format", choices=('export', 'eland_result'), default='export')
    optparser.add_option("-o", "--output",
                         metavar="FILE", help="write output to FILE (default output is to terminal)")
    
    # parse command line
    (options, args) = optparser.parse_args()
    # check options, args
    if len(args) < 1:
        optparser.error('No input file specified')
    
    # get arguments
    elandfilename = args[0]
    
    # setup output
    if options.output:
        outfile = open(options.output, 'w')
    else:
        outfile = sys.stdout


    '''
    treatments/controls: dict mapping sample names to flowcell/lane/fileformats

    key: string containing name of control
    value: list of tuples containing three fields -- flowcell, lanes, and fileformat
        flowcell   -- string e.g. '20E46AAXX'
        lanes      -- tuple e.g. (6,) or for paired-end, ('5_1', '5_2')
        fileformat -- currently accepts only 'export' 

    example:
        treatments = [('30351AAXX', ('2_1', '2_2'), 'export'),
                      ('30351AAXX', ('3_1', '3_2'), 'export')]

        controls = [('20E46AAXX', ('6',), 'export'),
                    ('20E6WAAXX', ('8',), 'export')]}
    '''
    # default options
    macs_options = {'--mfold': 10,
                    '--pvalue': '1e-4',
                    '--verbose': 3}
    hpeak_options = {'-fmin': 175,
                     '-fmax': 225,
                     '-sig': '1e-5',
                     '-wig': None,
                     '-ann': None}

    # create the app controllers    
    macs_app = MACS(options=macs_options)
    hpeak_app = HPeak(options=hpeak_options)

    # run
    for peakfinder in (macs_app,):
        find_peaks(peakfinder, vcap_erg, vcap_control)
        find_peaks(peakfinder, vcap_ku80, vcap_control)
        find_peaks(peakfinder, vcap_dnapk, vcap_control)
        find_peaks(peakfinder, vcap_ph2ax, vcap_control)
        find_peaks(peakfinder, vcap_foxa1, vcap_control)
        find_peaks(peakfinder, vcap_r1881_erg, vcap_control)